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Abstract. We study the effects of the exchange interaction between an adsorbed 
magnetic atom with easy-axis magnetic anisotropy and the conduction-band electrons 
from the substrate. We model the system using an anisotropic Kondo model and we 
compute the impurity spectral function which is related to the differential conductance 
(dl/dV) spectra measured using a scanning tunneling microscope. To make contact 
with the known experimental results for iron atoms on the CuN/Cu(100) surface 
[Hirjibehedin et al.. Science 317, 1199 (2007)], we calculated the spectral functions 
in the presence of an external magnetic field of varying strength applied along all three 
spatial directions. It is possible to establish an upper bound on the coupling constant 
J: in the range of the magnetic fields for which the experimental results are currently 
known (up to 7 T) , the low-energy features in the calculated spectra agree well with 
the measured dl/dV spectra if the exchange coupling constant J is at most half as 
large as that for cobalt atoms on the same surface. We show that for even higher 
magnetic field (between 8 and 9 T) applied along the "hollow direction" , the impurity 
energy states cross, giving rise to a Kondo effect which takes the form of a zero-bias 
resonance. The coupling strength J could be determined experimentally by performing 
tunneling spectroscopy in this range of magnetic fields. On the technical side, the paper 
introduces an approach for calculating the expectation values of global spin operators 
and all components of the impurity magnetic susceptibility tensor (including the out- 
of-diagonal ones) in numerical renormalization group (NRG) calculations with no spin 
symmetry. An appendix contains a dcnsity-functional-theory (DFT) study of the Co 
and Fe adsorbates on CuN/Cu(100) surface: we compare magnetic moments, as well 
as orbital energies, occupancies, centers, and spreads by calculating the maximally 
localized Wannier orbitals of the adsorbates. 



PACS numbers: 75.30.Gw, 72.10.Fk, 72.15.Qm 
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1. Introduction 

Recent advances in tlie experimental techniques, such as low-temperature scanning 
tunnehng microscopy (STM) and the X-ray magnetic circular dichroism (XMCD), made 
it possible to study magnetic properties of single adsorbed atoms on various surfaces 
[H El El m El El [7]. Particularly interesting are magnetic impurities on noble metal 
surfaces where many-particle effects such as the Kondo effect play an important role 
P|9]. Due to the reduced symmetry in the surface region, such magnetic adsorbates are 
strongly anisotropic [Il[10]. A simple parametrization of the leading magnetic anisotropy 
terms takes the form of 

H a.niso = DS^ + E{Sl — Sy), (1) 

where the S^, Sy, and are the cartesian components of the quantum- mechanical 
spin-S* operator, while the parameters D and E are known as the longitudinal and 
transverse magnetic anisotropy, respectively. If -D < (the "easy-axis" case), the spin 
tends to maximize the absolute value of the z component, thus the = ±S states 
dominate in the ground state doublet, which may be split for integer S even in the 
absence of an external magnetic field due to the transverse magnetic anisotropy E. If 
D > (the "easy-plane" or "hard-axis" case), the ground-state is a doublet consisting 
mostly of the 5*^ = ±1/2 states for half- integer spin, or a singlet consisting mostly of 
the Sz = state for integer spin. Similar models are also used to describe molecular 
magnets [III [121 [I3l [H [I5] . 

When adsorbed on metallic surfaces, the coupling of the spin to the conduction-band 
electrons can introduce new phenomena, the Kondo effect being the most interesting 
due to its unique signature expected in the tunneling spectra. For half-integer spins 
with easy-plane anisotropy, the level degeneracy in the ground state actually permits the 
occurrence of the Kondo effect at zero magnetic field, as indeed observed experimentally 
in the system of cobalt atoms adsorbed on the CuN ultra-thin layers deposited on the 
Cu(lOO) substrate surface [16]. For magnetic atoms with easy-axis anisotropy, however, 
the Kondo effect has not yet been observed. Its absence is either due to the zero-field 
level splitting caused by the transverse magnetic anisotropy E, which in the case of 
integer 5" plays the role of an effective magnetic field which suppresses the Kondo effect, 
or (for half- integer or for ii^ ^ \D\) due to the fact that the spin-flip scattering events 
of the conduction-band electrons can only change the impurity spin component by 1. 
Thus the two impurity ground-state levels only couple via higher-order processes and 
the Kondo screening is strongly suppressed, i.e., the Kondo scale is extremely small [T2] . 
Nevertheless, interesting many-particle effects in fact do occur in the case of easy-axis 
anisotropy, thus we study this class of the problems in the following parts of the paper, 
focusing in particular on the case of spin 2 as appropriate for iron atoms adsorbed on 
the CuN/Cu(100) substrate [10]. 

The absence of the Kondo effect at zero magnetic field does not imply, that 
the exchange coupling of an easy-axis magnetic impurity to the substrate is of no 
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consequence. In this work we show that the value of the exchange couphng constant J, 
defined through the Kondo couphng Hamiltonian 

Hk = Js- S, (2) 

where s is the spin-density of the conduction-band electrons at the position of the 
impurity, affects the impurity spectral function in a measurable way. The most drastic 
effect is the occurrence of a Kondo effect induced by a magnetic field, similar to the 
singlet-triplet Kondo effect in quantum dot systems |T7l HI [191 1201 HU HH |22] : for some 
high enough external magnetic field applied along the x-axis, the impurity level crossing 
leads to a situation where the conditions for the occurrence of a Kondo effect are met. 
We will show that the presently available experimental results (measured up to 7T, 
Ref. [in]) just stop short of the region where this Kondo resonance should be observed 
(between 8T and 9T). For this reason, the presently known experimental results can 
only be used to establish an upper bound on the value of J, but by extending the 
magnetic field range into the 9 T range, the Kondo resonance could be probed and the 
exchange coupling J might be determined even quantitatively. 

The paper is structured as follows. In Sec. [2] we provide the full definition of the 
model under study and describe the numerical technique which is used to compute the 
impurity spectral function. We also comment on the relation between the impurity 
spectral function and the experimentally measured dl/dV spectrum. In Sec. [3] we 
discuss the impurity level structure and the predicted dl/dV spectra in the limit of 
a decoupled impurity, J — )■ 0. Such spectra do not encompass any many-particle effects; 
however, they account accurately for the high-energy spin excitations which are observed 
experimentally. In Sec. H] the results of the numerical calculations for the full many- 
particle problem are presented and discussed in relation with the experimental results. 



Appendix A presents a comparative density-functional-theory (DFT) study of Co and 



Fe adatoms on the CuN/Cu(100) surface, while Appendix B contains the description of 
a method for computing the thermodynamic impurity spin susceptibility in problems 
having no symmetries in the spin space with the numerical renormalization group (NRG) 
method. 



2. Model and method 



We model the system using the Hamiltonian H = -ffband + Hk + -ffimp where 

Hha^nd = "^^kcla-Cka, (3) 
ka 

Hirap = gi^sB ■ S + -?/aniso- (4) 

The anisotropy term i^amso is given by Eq. ([1]), and the coupling term ifx by Eq. ([2]). 
The operator Cka- corresponds to a conduction-band electron with momentum k, spin 
a G {t, i}, and energy e^. In terms of these operators, the spin density of the conduction- 
band electrons at the position of the impurity (assumed to be at the origin of the 
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coordinate system) is given by s = (1/2)/q ^cra/3/0,/3, where 

foa = -yff^Cka, (5) 



is the number of the conduction-band states, and cr^^ is the vector of Pauh matrices. 
Finally, g is the gyromagnetic factor (assumed to be isotropic), and fis ~ 0.0579meV/T 
is the Bohr magneton. For an iron adatom on the CuN/Cu(100) surface, the following 
parameters have been established: spin S = 2, g = 2.11 ± 0.05, D = —1.55 ± 0.01 meV, 
E = 0.31 ± 0.01 meV (Ref. ^Q\)- For a complete description one should, furthermore, 
know the density of states (DOS) of the conduction band, p{u) = (l/N) — Cfc), 

and the value of the Kondo coupling J. In the interesting energy range (on the order 
of 10 meV around the Fermi level) the energy dependence of the DOS may be neglected 
and we can use a constant value p = p{u = 0). The low-energy behavior of the problem 
then depends solely on the dimensionless combination pj. We note that the model does 
not include any orbital moment degrees of freedom. This approximation is based on 
the fact that the orbital moment of an impurity adsorbed on the surface is typically 
strongly quenched due to the surface field effects (for the particular case of Fe and Co 
adsorbates on the surface of CuN/Cu(100), the total suppression of orbital degeneracy 
follows from the DFT results tabulated in Appendix A ). 

We solve the problem using the numerical renormalization group [231 1211 [25] (for the 
details on the application of the NRG to the problems of this class see also Refs. [261127] ). 
In this work, the focus of our interest will be the impurity spectral function defined as 

[2H] 

A^{u) = --lmT^{uj + i6), (6) 

TT 

where Tcr{u}) is the T-matrix for the scattering of electrons of spin a on the spin-S* 
impurity, given by 

Uuj) = ^{S,) + {{0^-Oi))^, (7) 
with the operators 

Ot = ^ {foiS- + /ot^,) , (8) 

= ^ {fotS^ - foiS.) . (9) 

Here = Sx ± iSy. The notation {{A;B))i^ denotes the Laplace transform of the 
correlator {{A;B))t = —i9{t){{A(t),B{0)}^) between two (fermionic) operators A and 
B, i.e., {{A;B))^ = e'^^^iA; B))tdt with Imw > 0; {A,B}+ = AB + BA is the 
ant i- commut at or . 

The T-matrix itself is defined through 

where Gkkai^) = ((cfco-; c|.^)) is the conduction-band electron Green's function and 
its unperturbed counterpart. In the STM experiments, the spin components are 
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not resolved, thus we will mostly plot the spin-averaged spectral functions, A{u) = 
[A^{u) + A^{uj)]/2. 

3. Decoupled impurity limit 

The spin excitation spectrum of a magnetic atom adsorbed on a thin insulating layer 
may for the most part be accounted for by fully neglecting its exchange coupling to the 
substrate electrons. According to this picture, the tunneling differential conductance 
will increase stepwise when the bias voltage is increased beyond the values {Ei — Egs)/e, 
where Ei are the energies of the excited states and Egg is the ground-state energy 
[21 [31 Ho] . At each of these voltages, an additional inelastic scattering (spin-flip) channel 
opens, typically resulting in an increased conductance |29l |30l [311 [321 |33l [Ml [351 [Ml [37] . 
In fact, the experiments have shown that the conductance step heights are described to 
a good approximation by p] 

mS,\gs)\' + mSy\gs)\'+mS,\gs)\^ (11) 

i.e., by the transition matrix elements for the spin operator. This empirical observation 
has recently received theoretical support [381 |39]. The expression suggests that the 
impurity is well thermalized with the substrate: the impurity relaxes to the ground-state 
on a time-scale which is short compared with the mean time between the tunneling 
events. Using this simple approach, one can compute dl/dV spectra which compare 
very well with the experimental results (as long as the many-particle effects are not 
important) if thermal broadening is taken into account. Here we show the results 
obtained by this procedure without adding any artificial thermal smearing, see Fig. [1^. 
The plots should be compared with Fig. 2A,B,F and Fig. S1,A,C in Ref. [lO]. Presenting 
the results in the zero-temperature limit uncovers more details, in particular the crossing 
of the the levels in the case of a field applied along the x-axis (this corresponds to the 
"hollow axis" in Ref. [10]). 

As is well known, the experimental results are in very good agreement with the 
simulated dl/dV spectra [10]. One comment, however, is in order. In experiments, the 
step height of the first step in the dl/dV spectra is observed to decrease as the field 
strength along the x-axis is increased. This can be explained, on one hand, by the 
thermal broadening of the spectra which leads to reduced heights of spectral features 
of the width smaller than the width of the broadening kernel. Another explanation is 
based on a possibility of a slight misalignment of the magnetic field from the direction 
of the actual x-axis of the magnetic atom. A simple calculation (Fig. [TJo) shows that 
even a small misalignment of 1 ° (which is quite likely to occur) can suppress the step 
height of the first excitation (at T = 0) without significantly modifying those at the 
other magnetic excitation energies. 

We now comment on the important observation that the energy difference between 
the two lowest levels decreases as the field strength along the x-axis is increased, see 
Fig. [It, top. A simple calculation neglecting any possible many-particle effects indicates 
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Figure 1. a) Simulated lETS spectra. No thermal broadening has been applied, b) 
Simulated lETS spectra for a field with a direction almost along the x-axis with a 
small component along the z-axis (misalignment by 1 ° . c) Energy-level diagrams for 
a field along the x-axis (top) and for a field misaligned by 1 ° (bottom). 



that the levels become degenerate at 

= -^J2EiE - D)^ 8.8 T. (12) 

For this field strength, a magnetic-field-induced Kondo effect can occur if the Kondo 
exchange coupling is antiferromagnetic, J > 0. We note that the aforementioned 
misalignment of the field will also tend to suppress the Kondo effect at ^ 8.8 T, 
since a component of the field along the z-axis will mix the levels, thereby leading to an 
avoided crossing of the levels (Fig. [It, bottom). In experiments seeking to detect the 
field-induced Kondo effect, it will be thus essential to not only fine-tune the magnitude 
of the magnetic field, but also its exact direction. 

4. Impurity spectral functions 

In Fig. [2] we plot the impurity spectral functions A{u) computed using the numerical 
renormalization group (NRG). The calculation is performed for zero temperature, thus 
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the widths of spectral features are due to intrinsic effects (the NRG spectral function 
overbroadening artifacts are small in these calculations; for discussion see Ref. |1Q]). 
Compared with the simulated dl/dV curves, all spectral features except for those near 
the Fermi level are almost completely smeared out. This is at odds with the experimental 
results, where the steps are observed to be sharp apart from the thermal broadening. 
While the calculated spectral function does encompass some inelastic scattering effects 
(see Fig. 3 in Ref. |27]), the tunneling current includes additional inelastic scattering 
contributions which are beyond the scope of our calculations that are constrained to 
the case of the magnetic impurity being in thermal equilibrium with the substrate 
electron reservoir at all times and do not include any additional interactions between the 
tunneling electron and the impurity local moment. Nevertheless, the results obtained 
from a NRG calculation provide a reliable description of the behavior of the system 
in the immediate vicinity of the Fermi level, which is probed by tunneling electrons 
of sufficiently low-energy so that the system is not driven significantly out of the 
equilibrium. 
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Figure 2. Impurity spectral functions A{uj) computed using the NRG. The 
horizontal energy axis has been transformed to show the equivalent bias voltages, 
while the vertical axis corresponds to the dl/dV differential conductance in the simple 
approximation where the LDOS is taken to be proportional to the spectral function. 
Here pJ = 0.05. 



The agreement between the low-energy features (interval [— ImV : +lmV]) in the 
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calculated and measured spectra PH] in the to 7T range is very good if thermal 
broadening is taken into account. The calculation was performed for pj = 0.05. This 
value is half as large as the exchange coupling for Co impurities on the same surface, 
pJ = 0.1, which has been established in Refs. [261 EZ]) by determining the value of 
pJ that allows to reproduce the experimental differential conductance spectra for the 
Co/CuN/Cu(100) adsorbate system. Due to similarities between Co and Fe atoms the 
exchange coupling constants should not be too different, but very likely of the same order 
of the magnitude and of the same sign (Co and Fe are neighbours in the periodic system; 
Appendix A for a more rigorous comparative study using the density functional 



see 



theory). In Fig. [3^ we compare the low-energy part of the spectral function for a fixed 
magnetic field of 7 T in the a;-axis direction for a range of values of J. For large J, 
the dip in the spectral density turns into a (Kondo) resonance which becomes fully 
developed for pJ > 0.1. For pJ in the interval 0.06 to 0.1, we find an intermediate 
structure: a dip with protrusions at the step edges. At high enough temperature, the 
protrusions will be washed out and a single antiresonance would be detected. Taking all 
these elements into account, the estimate of the upper bound pJ = 0.05 ±0.02 therefore 
seems reasonable. 
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Figure 3. Impurity spectral functions for 7 T magnetic field applied along the cc-axis 
for a range of pJ. 



The magnetic field at which the Kondo effect occurs for integer spin Fe is not given 
exactly by Eq. (fT2l) . but it rather depends on the value of J, since the different impurity 
levels are renormalized differently due to the anisotropy. In the J — )■ limit, the Kondo 
effect occurs at ~ 8.8 T, while for finite J the magnetic field where the Kondo peak 
may be observed is reduced. In addition, the parameter J controls the width of the 
Kondo resonance (and, of course, the Kondo temperature itself), thus for smaller J the 
interval of the magnetic fields where the Kondo effect may be observed is accordingly 
narrower and the Kondo scale lower. 

The magnetic-field-induced Kondo effect exhibits a number of rather peculiar 
features. We first note that it is very difficult to observe the fully developed Kondo 
resonance unless J is rather large. For the estimated exchange coupling of pJ = 0.05, 
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where the Kondo effect occurs for B* = 8.535 T, the Kondo temperature is actually 
extremely low, on the order of 0.1 mK, see Fig. HI Even at dilution refrigerator 
temperatures, it thus appears unlikely that a fully developed Kondo resonance will 
be observed. This does not, however, preclude the experimental detection of the 
Kondo physics at play in this system, since the evolution of the tunneling spectrum 
as the magnetic field is tuned should exhibit a characteristic merging of the two "half- 
resonances" . 

The results for the thermodynamic properties also reveal that the zz component 
of the magnetic susceptibility multiplied by the temperature [i.e., the effective moment 
Txzz{T)] exhibits a plateau at very high value of 3.5 (in suitable units), which is rather 
near the limiting value for a decoupled = ±2 doublet, which is 4. This plateau 
corresponds to a In 2 plateau in the impurity entropy, i.e., to an effective two-state 
system. On the scale of Tk, the moment is then screened. Note that the temperature 
dependence of TxzziT) is different (Kondo-like, see the arrow in Fig. H]) for = B* 
as compared with those for other values of the magnetic field, which merely follow the 
exp{—AE/kBT) rule, where AE is the difference between the two lowest impurity level 
energies. 

It is equally interesting to note that the spectral functions change discontinuously 
as we cross the transition value B*, see Fig. |5^; this corresponds to the system being 
in different ground states on either side of the transition. In the presence of a small 
magnetic field component along the z-axis, the level-crossing is replaced by a rapid 
cross-over. Any such rapid change of the spectral function over an extended energy 
interval would be an indicator of the proximity to the Kondo regime, even at higher 
temperatures. 

The magnetic-field-induced Kondo effect is not the only observable consequence 
of the exchange coupling with the substrate. The value (as well as the sign) of J is 
actually reflected in the spectral function for any value of the magnetic field. In Fig. [6] 
we compare, for example, the impurity spectral function at zero field for a range of 
coupling constants. The exchange coupling strongly affects the features around the 
Fermi level, but also those at higher energies, although to a lesser degree. 

The results which we have established for the S = 2 Kondo model also apply more 
generally for other integer spins: at some B* there will be level crossing which can induce 
a Kondo resonance at the Fermi level, but no such crossing occurs for the magnetic field 
along the y and z axis. The same is true for half-integer spin models, but in this case 
there is, in addition, a two-fold level degeneracy also for zero magnetic field, thus for 
half-integer-spin magnetic impurities one could observe both the zero-field Kondo effect 
and the magnetic-field- induced Kondo effect. 

5. Conclusion 

We have studied thermodynamic and spectral properties of the Kondo model with an 
easy-axis anisotropy. By comparing the calculated spectral functions with the measured 
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kg T (meV) 

Figure 4. Thermodynamic properties for a range of magnetic fields along the x-axis. 
We plot the impurity contribution to the magnetic susceptibility in various directions, 
Ximp,aa, where a € {x,y,z} is a direction in space, and the impurity contribution to 
the entropy, S'imp- The Kondo temperature is estimated from the data for S'imp to be 
Tk ~ 10~^ meV ^ 0.1 mK. The arrow points to the curve with exhibits Kondo-like 
temperature dependence. 

ones in the magnetic field range up to 7 T we have estimated the upper bound for the 
exchange couphng of a Fe adatom on CuN/Cu(100) surface with the bulk conduction- 
band electrons to be pj = 0.05 ± 0.02. By extending the range of the magnetic fields in 
the experiments up to 9 T, it should be possible to obtain an improved upper bound on 
J or perhaps even observe a magnetic-field-induced Kondo resonance if J is high enough. 
We note here that, strictly speaking, the quantity J scales under the renormalizat ion- 
group flow, thus its actual value depends on the choice of the bandwidth and the 
corresponding density of states p; a more rigorous treatment of the problem would take 
into account the full density-of-states dependence of the substrate and, furthermore, 
the momentum dependence of the Kondo coupling J{k,k'), see for example Ref |¥T] . 
Nevertheless, a simple characterization in terms of the dimensionless combination pJ is 
sufficient to compare the over-all effect of the Kondo screening for different adsorbed 
atoms. 

For quantitative estimates, it will be important also to reduce the temperature at 
which the experiments are performed as far as possible. Irrespective if the Kondo 
resonance around 8.5 T will eventually be observed or not, the real motivation for 
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Figure 5. Impurity spectral functions in the vicinity of the magnetic- field-induced 
Kendo effect at = B*. a) Overview plot, b) close-up on the low- voltage region. 
The jaggedness of the curve at the transition point is due to increased numeric artifacts 
related to the ground-state crossing, c) Spin-resolved spectral function away from the 
transition point. Note the rather sharp features in the field-split Kondo resonance. 




Figure 6. Impurity spectral functions at zero field for different values of exchange 
coupling J. 
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performing such experiments actually comes from the need to estimate the exchange 
coupling to the bulk electrons, since this parameter also determines the spin relaxation 
time; thus it plays an important role in the possible applications of high-spin states of 
magnetic impurities for storing and processing information. 

Appendix A. DFT study of Co and Fe adatoms on CuN/Cu(100) 



Table Al. Properties of the Co and Fe adatom on CuN/Cu(100) surface. 



System 


Moment (/is /cell) 


nd 




ris 


Co/Cu2N/Cu(100) 


2.88 


7.56 


2.10 


0.64 


Fe/Cu2N/Cu(100) 


4.12 


6.61 


3.09 


0.63 



Table A2. The properties of the MLWFs of Co on CuN/Cu(100). We tabulate orbital 
energy e, occupancy n, orbital spread s, and orbital center along the z-axis. The orbital 
labels correspond to the initial projections used in the calculation of the corresponding 
maximally localized Wannier orbitals, although the orbitals are significantly modified 
from the free-atom orbitals due to the strong hybridization with the neighboring N 
atoms and the Cu atom just below the adsorbate. The position of the Co adatom in 
the z-direction is 2;(Co) =-0.148 A; the origin for the z-axis corresponds to the top-most 
layer of Cu atoms prior to relaxation. 



Co Minority spin \. Majority spin t 



Orbital 


e (eV) 


n 




z(A) 


e(eV) 


n 


.(A^) 


z(A) 


s 


1.70 


0.27 


2.36 


-3.06 


1.32 


0.37 


2.07 


-3.03 




-1.01 


0.88 


0.73 


-2.01 


-2.41 


0.97 


0.52 


-2.06 


dzx 


-0.07 


0.28 


0.94 


-2.09 


-2.83 


1.00 


0.42 


-2.15 


dzy 


-0.50 


0.41 


0.50 


-2.17 


-2.98 


0.97 


0.44 


-2.17 




-1.04 


0.80 


1.16 


-2.03 


-2.30 


0.90 


0.77 


-2.08 


dxy 


-0.34 


0.36 


0.65 


-2.12 


-2.90 


0.99 


0.46 


-2.13 


d total 




2.73 








4.83 







To further substantiate our claim of the similarities between Co and Fe, 
which are neighbours in the periodic system differing by a single electron, we 
have performed a detailed comparative study using the density functional theory 
(DFT) [l2l US]. We used the PWSCF code (from the Quantum Espresso package. 



http://www.quantum-espresso.org/, Ref. [H]) with a plane-wave basis set, ultra-soft 



pseudopotentials, and Perdew-Burke-Ernzerhof cx-GGA exchange-correlation functional 
[ISl HHl m] . The kinetic-energy cutoff was 35 Ry and the density cut-off 400 Ry. We used 
an 8x8x4 Monkhorst-Pack mesh of /c-points in the self-consistent-field band-structure 
calculation (but 6x6x1 for the initial relaxation calculation). The cold-smearing by 
0.035 Ry has been performed using the Marzari-Vanderbilt scheme [17]. We used the 
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Table A3. The properties of the MLWFs of Fe on CuN/Cu(100). The position of the 
Fe adatom in the z-direction is 2(Fe) =-0.158 A. 



Fe 



Minority spin 4- 



Majority spin f 



Orbital e (eV) 



n 



s (A2) z (A) e (eV) 



n 



S (A2) Z (A) 




S 



1.59 
0.31 
0.63 
0.22 
0.30 
0.52 



0.22 2.34 

0.49 0.76 

0.10 0.98 

0.26 0.50 

0.72 1.29 

0.19 0.47 
1.76 



3.13 
2.18 
2.12 
2.35 
2.17 
2.27 



0.50 
3.26 
3.56 
3.33 
2.65 
3.19 



0.41 2.02 

0.97 0.47 

1.00 0.43 

0.99 0.51 

0.90 0.70 

0.99 0.47 
4.85 



3.16 
2.23 
2.30 
2.33 
2.23 
2.27 



experimentally-determined lattice constant for Cu, a = 0.361 nm. The substrate was 
modeled using a 2-by-2 supercell in the lateral direction consisting of 4 Cu layers (i.e., 37 
atoms in total including the N atoms and the adsorbate). The Cu atoms in the bottom- 
most layer were fixed, while all others were allowed to relax. We note that similar 
calculations have already been performed for Fe and Mn adatoms [10] and recently also 
for Mn chains |18] on the same surface. 

To gain more insight into the orbital structure of the impurities, we have also 
calculated the maximally localized Wannier functions (MLWF) |49l |50] on the impurity 
site using the WannierQO package (http:/ /ww w .wannier. org/ ^ Ref. [51]). The MLWFs 
provide a very compact and accurate local representation of the electronic structure: 
they give direct insight into the nature of the chemical bonding. A 4x4x2 uniform grid 
of fc-points was used here. The initial projections consisted of d and s states on the 
adatom and the Cu atoms, and of sp3 hybridized states on the N atom; in addition, 
s states were added in the interstitial sites of the lattice to describe the delocalized 
electrons (1 per Cu atom). The energy window was chosen so as to encompass all d and 
s levels of the adatom. 

Table lAll presents an overview of the results: the magnetic moment and the 
occupancies of the adatom levels. The results are expected - the two adatoms differ by 
one electron in the d orbital. More details can be found in Tables IA2I and I A3 1 where 
we show the orbital energies (defined as the average energy of the spectral function of 
a particular MLWF), occupancies (defined as the integral of the spectral function of 
MLWF up to the Fermi energy), the centers of the MLWFs, as well as their spreads. It 
is manifest that the difference consists of one additional electron in the "minority" spin 
orbitals for the Co adatom. In fact, the spectral functions shown in Fig. lAl I reveal that 
the "majority" spin spectral functions are essentially the same in both cases, while the 
"minority" spin spectral functions are to a good approximation just rigidly shifted. 
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Figure Al. Partial density of states (spectral functions) for the MLWF on Co (left) 
and Fe (right) adatoms. We plot both spin projections: "majority" spin (above), 
"minority" spin (below). The spectral functions were computed on a 16x16x3 grid of 
fc-points with a 0.1 eV binning resolution. 
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Appendix B. Calculation of the magnetic susceptibility in anisotropic 
models 

In general, the magnetic susceptibility is calculated as [521 [25] : 

Xa{T) = {gfiB? [ dr ((5,(r)5„(0)) - {S^f) , (B.l) 

where Sa is the a G {x,y,z} component of the spin and (3 = l/ksT. The impurity 
contribution to the total magnetic susceptibility, Ximp,«, is then obtained by subtracting 
the susceptibility of a system without the impurity, XaK'^)^ from the magnetic 
susceptibility of the full system (i.e., the continuum of the conductance-band electrons 
plus the impurity), XfuwAT)- 

XimpAT) = XfuuAT) - xfiT). (B.2) 

In the problems with U(l) spin symmetry, Sz is a conserved quantum number and 
it is used to classify the invariant subspaces of the Hilbert space, thus the calculation 
of the expectation values in fIB.ll) is particularly simple for a = z: 

X.{T) = {gf^EfP {{Sim - [{S.mf) , (B.3) 

where (5',)(T) = (1/Z(^)) Es,,™ exp(-/3E£^^) and similarly for (52)(T). Here Z(^) is 
the partition function at shell N, defined as Z^^^ = Tr7v(e"^^) = J2s^,n^w{~I^Eg^}^), 
Sz the conserved spin quantum number, n a multi-index composed of all the remaining 
quantum numbers used to classify the states, and E^^}^ the energy of the state with given 
quantum numbers. For other directions a, and for more general problems without any 
symmetry in the spin space, a full dynamical calculation of the correlator {Sa{T)Sa{0)) 
has to be performed. 

The calculation of Xa requires an iteration of the operator matrix elements for all 
components of the total spin operator Sa- The approach is similar to the calculation of 
the matrix elements of an arbitrary local impurity operator, but it is extended to allow 
for total ("global") spin operators Sa- The matrix elements at a given step + 1 are 
calculated from those at the previous step by using the unitary basis transformation 
resulting from the diagonalization of the Hamiltonian matrices. We follow the notation 
of Ref. [25] ; see Eqs. (39,40,43) in the cited work: 

\w)n+i = '^U{w,rs)\r;s)N+i, (B.4) 

rs 

with 

\r;s)N+i = \r)N®\siN+l)). (B.5) 

Here r indexes the eigenpairs of the chain Hamiltonian at the previous iteration, i.e., 
HN\r)N = Ef^{r)\r)N, the states |s(A^ + 1)) constitute a basis for the added site, 
while U{w, rs) is the unitary matrix which diagonalizes the Hamiltonian i^Ar+i, so that 

Hn^i\w) N+i = En+iIw) N+1- We write Sa~^ = Sa + s^"^^, where S is the total spin 
operator up to (and including) the site N, while s^~^^ denotes the local spin operator 
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on site + 1. The matrix elements are then updated as follows [compare with (75) in 
Ref. gg]: 

= n+i{w\S^+'\w')n+i (B.6) 
= E U{w, rs);v+i(r; s\S^+' ■ E U{w', r's') \r'; s^+i (B.7) 

rs r's' 

= J2 U{w,rs)U{w',r's'){M{r\S^\r)NSss' (B.8) 

rs,r's' 

+ 6rAs\s^^'\s')) (B.9) 
= T.T.U{w,rs)M^/U{w,r's) (B.IO) 

+ E f E U{w, rs)U{w', rs')) {s\s^\s'). (B.ll) 

ss' \ r / 

The first part of the expression on the right-hand side, (IB.IOI) . is exactly the same as 
the conventional prescription for the recursive calculation of the matrix elements of the 
operators which are singlets with respect to the symmetry group of the Hamiltonian: 
for each contribution s, one performs a similarity transformation of the matrix M^'^ 
from the previous iteration. The second part of the expression, ( IB. Ill) , consists in 
adding the matrix elements of the spin operator for the additional site weighted by the 
inner-products between the columns of submatrices of the unitary transformation U. 

At each iteration, the required expectation values are calculated in the same way 
as for local singlet operators [cf. (58) in Ref. [25]]: 

(5a)(T^)^^E^l (B.12) 

where the matrix elements M^^ are those that correspond to the operator S^- 

The imaginary-time integral of the correlator C(r) = {Sa{T)Sa{0)) may be 
expressed using the Lehmann representation as 

Note that for the spin operator along the y direction, the NRG calculation has to be 
performed using complex numbers. 

Where comparisons can be made (e.g., in isotropic models) the magnetic 
susceptibility calculated in the traditional way and using the proposed technique differ 
by no more than a few per mil in calculations performed for A = 2, and even less for 
higher A. We observed, however, that more states need to be kept in this approach 
than in the standard NRG thermodynamic calculations. The calculation of x is made 
possible by the property of the energy-scale separation of total spin operators, much 
like the entire NRG iteration is based on the energy-scale separation property of the 
Hamiltonian. 

The approach described in this section can be immediately generalized to the 
calculation of the out-of-diagonal components of the magnetic susceptibility tensor. 
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